From 2a64d98d0211d9045fe342f703294648f6760a1f Mon Sep 17 00:00:00 2001
From: Kutalmis Bercin <kutalmis.bercin@esi-group.com>
Date: Fri, 14 Aug 2020 12:50:15 +0100
Subject: [PATCH] ENH: PatchParticleHistogram: add a new cloud FO

  Computes a histogram for the distribution of particle diameters
  and corresponding number of particles hitting on a given list of patches.

  A minimal example by using `constant/reactingCloud1Properties.cloudFunctions`:

  ```
  patchParticleHistogram1
  {
      // Mandatory entries (unmodifiable)
      type             patchParticleHistogram;
      patches          (<patch1> <patch2> ... <patchN>);
      nBins            10;
      min              0.1;
      max              10.0;
      maxStoredParcels 20;
  }
  ```
---
 .../include/makeParcelCloudFunctionObjects.H  |   3 +
 .../makeReactingParcelCloudFunctionObjects.H  |   2 +
 .../PatchParticleHistogram.C                  | 254 ++++++++++++++++++
 .../PatchParticleHistogram.H                  | 222 +++++++++++++++
 .../filter/constant/reactingCloud1Properties  |  14 +
 5 files changed, 495 insertions(+)
 create mode 100644 src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchParticleHistogram/PatchParticleHistogram.C
 create mode 100644 src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchParticleHistogram/PatchParticleHistogram.H

diff --git a/src/lagrangian/intermediate/parcels/include/makeParcelCloudFunctionObjects.H b/src/lagrangian/intermediate/parcels/include/makeParcelCloudFunctionObjects.H
index 4615278aa0..ed6c4f6bf9 100644
--- a/src/lagrangian/intermediate/parcels/include/makeParcelCloudFunctionObjects.H
+++ b/src/lagrangian/intermediate/parcels/include/makeParcelCloudFunctionObjects.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2018 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,6 +38,7 @@ License
 #include "ParticleTrap.H"
 #include "PatchCollisionDensity.H"
 #include "PatchPostProcessing.H"
+#include "PatchParticleHistogram.H"
 #include "RemoveParcels.H"
 #include "VoidFraction.H"
 
@@ -53,6 +55,7 @@ License
     makeCloudFunctionObjectType(ParticleTrap, CloudType);                      \
     makeCloudFunctionObjectType(PatchCollisionDensity, CloudType);             \
     makeCloudFunctionObjectType(PatchPostProcessing, CloudType);               \
+    makeCloudFunctionObjectType(PatchParticleHistogram, CloudType);            \
     makeCloudFunctionObjectType(RemoveParcels, CloudType);                     \
     makeCloudFunctionObjectType(VoidFraction, CloudType);
 
diff --git a/src/lagrangian/intermediate/parcels/include/makeReactingParcelCloudFunctionObjects.H b/src/lagrangian/intermediate/parcels/include/makeReactingParcelCloudFunctionObjects.H
index a11f3de2e2..dd45507c68 100644
--- a/src/lagrangian/intermediate/parcels/include/makeReactingParcelCloudFunctionObjects.H
+++ b/src/lagrangian/intermediate/parcels/include/makeReactingParcelCloudFunctionObjects.H
@@ -38,6 +38,7 @@ License
 #include "ParticleTrap.H"
 #include "PatchCollisionDensity.H"
 #include "PatchPostProcessing.H"
+#include "PatchParticleHistogram.H"
 #include "RemoveParcels.H"
 #include "VoidFraction.H"
 #include "WeberNumberReacting.H"
@@ -55,6 +56,7 @@ License
     makeCloudFunctionObjectType(ParticleTrap, CloudType);                      \
     makeCloudFunctionObjectType(PatchCollisionDensity, CloudType);             \
     makeCloudFunctionObjectType(PatchPostProcessing, CloudType);               \
+    makeCloudFunctionObjectType(PatchParticleHistogram, CloudType);            \
     makeCloudFunctionObjectType(RemoveParcels, CloudType);                     \
     makeCloudFunctionObjectType(VoidFraction, CloudType);                      \
     makeCloudFunctionObjectType(WeberNumberReacting, CloudType);
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchParticleHistogram/PatchParticleHistogram.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchParticleHistogram/PatchParticleHistogram.C
new file mode 100644
index 0000000000..d60f1723dc
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchParticleHistogram/PatchParticleHistogram.C
@@ -0,0 +1,254 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "PatchParticleHistogram.H"
+#include "Pstream.H"
+#include "stringListOps.H"
+#include "ListOps.H"
+#include "ListListOps.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::label Foam::PatchParticleHistogram<CloudType>::applyToPatch
+(
+    const label globalPatchi
+) const
+{
+    return patchIDs_.find(globalPatchi);
+}
+
+
+// * * * * * * * * * * * * * protected Member Functions  * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::PatchParticleHistogram<CloudType>::write()
+{
+    forAll(times_, i)
+    {
+        List<List<scalar>> procTimes(Pstream::nProcs());
+        procTimes[Pstream::myProcNo()] = times_[i];
+        Pstream::gatherList(procTimes);
+
+        List<List<scalar>> procDiameters(Pstream::nProcs());
+        procDiameters[Pstream::myProcNo()] = patchDiameters_[i];
+        Pstream::gatherList(procDiameters);
+
+        List<List<scalar>> procParticles(Pstream::nProcs());
+        procParticles[Pstream::myProcNo()] = patchParticles_[i];
+        Pstream::gatherList(procParticles);
+
+        if (Pstream::master())
+        {
+            const fvMesh& mesh = this->owner().mesh();
+
+            mkDir(this->writeTimeDir());
+
+            const word& patchName = mesh.boundaryMesh()[patchIDs_[i]].name();
+
+            OFstream patchOutFile
+            (
+                this->writeTimeDir()/patchName + ".post",
+                IOstream::ASCII,
+                IOstream::currentVersion,
+                mesh.time().writeCompression()
+            );
+
+            List<scalar> globalTimes;
+            globalTimes = ListListOps::combine<List<scalar>>
+            (
+                procTimes,
+                accessOp<List<scalar>>()
+            );
+
+            List<scalar> globalDiameters;
+            globalDiameters = ListListOps::combine<List<scalar>>
+            (
+                procDiameters,
+                accessOp<List<scalar>>()
+            );
+
+            List<scalar> globalParticles;
+            globalParticles = ListListOps::combine<List<scalar>>
+            (
+                procParticles,
+                accessOp<List<scalar>>()
+            );
+
+            // Compute histogram
+            List<scalar> nParticles(nBins_, Zero);
+            forAll(globalDiameters, j)
+            {
+                const label bini = (globalDiameters[j] - min_)/delta_;
+                if (bini >= 0 && bini < nBins_)
+                {
+                    nParticles[bini] += globalParticles[j];
+                    nParticlesCumulative_[i][bini] += globalParticles[j];
+                }
+            }
+
+            patchOutFile
+                << "# nBin=" << nBins_
+                << "; min="  << min_
+                << "; max="  << max_ << nl
+                << "# d            nParticles            nParticlesCumulative"
+                << endl;
+
+            forAll(nParticles, j)
+            {
+                patchOutFile
+                    << binEdges_[j]
+                    << "-"
+                    << binEdges_[j + 1]
+                    << " "
+                    << nParticles[j]
+                    << " "
+                    << nParticlesCumulative_[i][j]
+                    << nl;
+            }
+        }
+
+        times_[i].clearStorage();
+        patchDiameters_[i].clearStorage();
+        patchParticles_[i].clearStorage();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::PatchParticleHistogram<CloudType>::PatchParticleHistogram
+(
+    const dictionary& dict,
+    CloudType& owner,
+    const word& modelName
+)
+:
+    CloudFunctionObject<CloudType>(dict, owner, modelName, typeName),
+    nBins_(dict.getCheck<label>("nBins", labelMinMax::ge(1))),
+    min_(dict.getScalar("min")),
+    max_(dict.getScalar("max")),
+    delta_((max_ - min_)/scalar(nBins_)),
+    maxStoredParcels_(dict.getScalar("maxStoredParcels")),
+    binEdges_(nBins_ + 1),
+    patchIDs_(),
+    nParticlesCumulative_(),
+    times_(),
+    patchDiameters_(),
+    patchParticles_()
+{
+    if (min_ >= max_)
+    {
+        FatalIOErrorInFunction(dict)
+            << "Histogram minimum = " << min_
+            << ", cannot be larger than histogram maximum = " << max_
+            << exit(FatalIOError);
+    }
+
+    if (maxStoredParcels_ <= 0)
+    {
+        FatalIOErrorInFunction(dict)
+            << "maxStoredParcels = " << maxStoredParcels_
+            << ", cannot be equal to or less than zero"
+            << exit(FatalIOError);
+    }
+
+    // Compute histogram-bin properties
+    binEdges_[0] = min_;
+    for (label i = 0; i < nBins_; ++i)
+    {
+        const scalar next = min_ + (i+1)*delta_;
+        binEdges_[i+1] = next;
+    }
+
+    // Compute histogram-patch properties
+    const wordRes patchMatcher(dict.get<wordRes>("patches"));
+
+    patchIDs_ = patchMatcher.matching(owner.mesh().boundaryMesh().names());
+
+    if (patchIDs_.empty())
+    {
+        FatalIOErrorInFunction(dict)
+            << "No matching patches found: "
+            << flatOutput(patchMatcher) << nl
+            << exit(FatalIOError);
+    }
+
+    nParticlesCumulative_ =
+        List<List<scalar>>(patchIDs_.size(), List<scalar>(nBins_, Zero));
+
+    times_.setSize(patchIDs_.size());
+    patchDiameters_.setSize(patchIDs_.size());
+    patchParticles_.setSize(patchIDs_.size());
+}
+
+
+template<class CloudType>
+Foam::PatchParticleHistogram<CloudType>::PatchParticleHistogram
+(
+    const PatchParticleHistogram<CloudType>& ppm
+)
+:
+    CloudFunctionObject<CloudType>(ppm),
+    nBins_(ppm.nBins_),
+    min_(ppm.min_),
+    max_(ppm.max_),
+    delta_(ppm.delta_),
+    maxStoredParcels_(ppm.maxStoredParcels_),
+    binEdges_(ppm.binEdges_),
+    patchIDs_(ppm.patchIDs_),
+    nParticlesCumulative_(ppm.nParticlesCumulative_),
+    times_(ppm.times_),
+    patchDiameters_(ppm.patchDiameters_),
+    patchParticles_(ppm.patchParticles_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::PatchParticleHistogram<CloudType>::postPatch
+(
+    const parcelType& p,
+    const polyPatch& pp,
+    bool&
+)
+{
+    const label patchi = pp.index();
+    const label localPatchi = applyToPatch(patchi);
+
+    if (localPatchi != -1 && times_[localPatchi].size() < maxStoredParcels_)
+    {
+        times_[localPatchi].append(this->owner().time().value());
+        patchDiameters_[localPatchi].append(p.d());
+        patchParticles_[localPatchi].append(p.nParticle());
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchParticleHistogram/PatchParticleHistogram.H b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchParticleHistogram/PatchParticleHistogram.H
new file mode 100644
index 0000000000..c8b3b99239
--- /dev/null
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchParticleHistogram/PatchParticleHistogram.H
@@ -0,0 +1,222 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::PatchParticleHistogram
+
+Description
+    Computes a histogram for the distribution of particle diameters
+    and corresponding number of particles hitting on a given list of patches.
+
+    Operands:
+    \table
+      Operand        | Type | Location
+      input          | -    | -
+      output file    | dat  | $FOAM_CASE/postProcessing/\<FO\>/\<time\>/\<file\>
+      output field   | -    | -
+    \endtable
+
+    The output file contains two columns, the first is the bin edges of
+    the particle diameter (i.e. \c d), and the second is the number of
+    particles whose diameter falling into the corresponding bin
+    (i.e. \c nParticles).
+
+Usage
+    Minimal example by using
+    \c constant/reactingCloud1Properties.cloudFunctions:
+    \verbatim
+    patchParticleHistogram1
+    {
+        // Mandatory entries (unmodifiable)
+        type             patchParticleHistogram;
+        patches          (<patch1> <patch2> ... <patchN>);
+        nBins            10;
+        min              0.1;
+        max              10.0;
+        maxStoredParcels 20;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                        | Type   | Reqd  | Dflt
+      type         | Type name: patchParticleHistogram  | word   |  yes  | -
+      patches      | Names of operand patches         | wordList |  yes  | -
+      nBins        | Number of histogram bins           | label  |  yes  | -
+      max          | Maximum value of histogram data    | scalar |  yes  | -
+      min          | Minimum value of histogram data    | scalar |  yes  | -
+      maxStoredParcels | Maximum number of parcels to process | label | yes | -
+    \endtable
+
+Note
+  - The underlying type of \c maxStoredParcels is set as a scalar for I/O.
+
+SourceFiles
+    PatchParticleHistogram.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PatchParticleHistogram_H
+#define PatchParticleHistogram_H
+
+#include "CloudFunctionObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class PatchParticleHistogram Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class PatchParticleHistogram
+:
+    public CloudFunctionObject<CloudType>
+{
+    // Private Data
+
+        //- Convenience typedef for parcel type
+        typedef typename CloudType::particleType parcelType;
+
+        //- Number of data bins
+        const label nBins_;
+
+        //- Minimum value of histogram data
+        const scalar min_;
+
+        //- Maximum value of histogram data
+        const scalar max_;
+
+        //- Bin width of histogram
+        const scalar delta_;
+
+        //- Maximum number of parcels to store - set as a scalar for I/O
+        const scalar maxStoredParcels_;
+
+        //- Bin edges of histogram
+        scalarField binEdges_;
+
+        //- List of patch indices to post-process
+        labelList patchIDs_;
+
+        //- Accumulated number of particles per patch
+        //- binned according to the histogram settings
+        List<List<scalar>> nParticlesCumulative_;
+
+        //- List of time for each data record
+        List<DynamicList<scalar>> times_;
+
+        // List of patch-hit particle diameters
+        List<DynamicList<scalar>> patchDiameters_;
+
+        // List of number of patch-hit particles
+        List<DynamicList<scalar>> patchParticles_;
+
+
+    // Private Member Functions
+
+        //- Return local patchi if patch is in patchIds_ list
+        label applyToPatch(const label globalPatchi) const;
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Write post-processing info
+        void write();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("patchParticleHistogram");
+
+
+    // Constructors
+
+        //- No default construct
+        PatchParticleHistogram() = delete;
+
+        //- Construct from dictionary
+        PatchParticleHistogram
+        (
+            const dictionary& dict,
+            CloudType& owner,
+            const word& modelName
+        );
+
+        //- Copy construct
+        PatchParticleHistogram(const PatchParticleHistogram<CloudType>& ppm);
+
+        //- No copy assignment
+        void operator=(const PatchParticleHistogram<CloudType>&) = delete;
+
+        //- Construct and return a clone
+        virtual autoPtr<CloudFunctionObject<CloudType>> clone() const
+        {
+            return autoPtr<CloudFunctionObject<CloudType>>
+            (
+                new PatchParticleHistogram<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~PatchParticleHistogram() = default;
+
+
+    // Member Functions
+
+        // Evaluation
+
+            //- Post-patch hook
+            virtual void postPatch
+            (
+                const parcelType& p,
+                const polyPatch& pp,
+                bool& keepParticle
+            );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "PatchParticleHistogram.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFoam/filter/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingParcelFoam/filter/constant/reactingCloud1Properties
index 7b93408646..bbda9f138c 100644
--- a/tutorials/lagrangian/reactingParcelFoam/filter/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/reactingParcelFoam/filter/constant/reactingCloud1Properties
@@ -161,6 +161,20 @@ subModels
 
 cloudFunctions
 {
+    patchParticleHistogram1
+    {
+        type            patchParticleHistogram;
+        patches
+        (
+            cycLeft_half0
+            cycLeft_half1
+        );
+        nBins            30;
+        min              0.0009;
+        max              0.001;
+        maxStoredParcels 20;
+    }
+
     patchPostProcessing1
     {
         type            patchPostProcessing;