From 611959fae118edb259fe38bf766d24c758c74723 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Fri, 12 Jul 2013 13:14:55 +0100
Subject: [PATCH] ENH: PatchflowRateInjection model - refactored and improved

---
 .../PatchFlowRateInjection.C                  | 192 ++++++------------
 .../PatchFlowRateInjection.H                  |  49 ++---
 2 files changed, 81 insertions(+), 160 deletions(-)

diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchFlowRateInjection/PatchFlowRateInjection.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchFlowRateInjection/PatchFlowRateInjection.C
index af9fd37d03..8499589784 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchFlowRateInjection/PatchFlowRateInjection.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchFlowRateInjection/PatchFlowRateInjection.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,6 +27,7 @@ License
 #include "TimeDataEntry.H"
 #include "distributionModel.H"
 #include "mathematicalConstants.H"
+#include "surfaceFields.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -39,19 +40,23 @@ Foam::PatchFlowRateInjection<CloudType>::PatchFlowRateInjection
 )
 :
     InjectionModel<CloudType>(dict, owner, modelName,typeName),
-    patchName_(this->coeffDict().lookup("patchName")),
-    patchId_(owner.mesh().boundaryMesh().findPatchID(patchName_)),
-    patchArea_(0.0),
-    patchNormal_(vector::zero),
+    patchInjectionBase(owner.mesh(), this->coeffDict().lookup("patchName")),
     phiName_(this->coeffDict().template lookupOrDefault<word>("phi", "phi")),
     rhoName_(this->coeffDict().template lookupOrDefault<word>("rho", "rho")),
     duration_(readScalar(this->coeffDict().lookup("duration"))),
-    concentration_(readScalar(this->coeffDict().lookup("concentration"))),
-    parcelsPerSecond_
+    concentration_
     (
-        readScalar(this->coeffDict().lookup("parcelsPerSecond"))
+        TimeDataEntry<scalar>
+        (
+            owner.db().time(),
+            "concentration",
+            this->coeffDict()
+        )
+    ),
+    parcelConcentration_
+    (
+        readScalar(this->coeffDict().lookup("parcelConcentration"))
     ),
-    U0_(vector::zero),
     sizeDistribution_
     (
         distributionModels::distributionModel::New
@@ -59,45 +64,11 @@ Foam::PatchFlowRateInjection<CloudType>::PatchFlowRateInjection
             this->coeffDict().subDict("sizeDistribution"),
             owner.rndGen()
         )
-    ),
-    cellOwners_(),
-    fraction_(1.0),
-    pMeanVolume_(0.0)
+    )
 {
-    if (patchId_ < 0)
-    {
-        FatalErrorIn
-        (
-            "PatchFlowRateInjection<CloudType>::PatchFlowRateInjection"
-            "("
-                "const dictionary&, "
-                "CloudType&"
-            ")"
-        )   << "Requested patch " << patchName_ << " not found" << nl
-            << "Available patches are: " << owner.mesh().boundaryMesh().names()
-            << nl << exit(FatalError);
-    }
-
-    const polyPatch& patch = owner.mesh().boundaryMesh()[patchId_];
-
     duration_ = owner.db().time().userTimeToTime(duration_);
 
-    updateMesh();
-
-    // TODO: retrieve mean diameter from distrution model
-    scalar pMeanDiameter =
-        readScalar(this->coeffDict().lookup("meanParticleDiameter"));
-    pMeanVolume_ = constant::mathematical::pi*pow3(pMeanDiameter)/6.0;
-
-    // patch geometry
-    label patchSize = cellOwners_.size();
-    label totalPatchSize = patchSize;
-    reduce(totalPatchSize, sumOp<label>());
-    fraction_ = scalar(patchSize)/totalPatchSize;
-
-    patchArea_ = gSum(mag(patch.faceAreas()));
-    patchNormal_ = gSum(patch.faceNormals())/totalPatchSize;
-    patchNormal_ /= mag(patchNormal_);
+    patchInjectionBase::updateMesh(owner.mesh());
 
     // Re-initialise total mass/volume to inject to zero
     // - will be reset during each injection
@@ -113,20 +84,13 @@ Foam::PatchFlowRateInjection<CloudType>::PatchFlowRateInjection
 )
 :
     InjectionModel<CloudType>(im),
-    patchName_(im.patchName_),
-    patchId_(im.patchId_),
-    patchArea_(im.patchArea_),
-    patchNormal_(im.patchNormal_),
+    patchInjectionBase(im),
     phiName_(im.phiName_),
     rhoName_(im.rhoName_),
     duration_(im.duration_),
     concentration_(im.concentration_),
-    parcelsPerSecond_(im.parcelsPerSecond_),
-    U0_(im.U0_),
-    sizeDistribution_(im.sizeDistribution_().clone().ptr()),
-    cellOwners_(im.cellOwners_),
-    fraction_(im.fraction_),
-    pMeanVolume_(im.pMeanVolume_)
+    parcelConcentration_(im.parcelConcentration_),
+    sizeDistribution_(im.sizeDistribution_().clone().ptr())
 {}
 
 
@@ -142,9 +106,7 @@ Foam::PatchFlowRateInjection<CloudType>::~PatchFlowRateInjection()
 template<class CloudType>
 void Foam::PatchFlowRateInjection<CloudType>::updateMesh()
 {
-    // Set/cache the injector cells
-    const polyPatch& patch = this->owner().mesh().boundaryMesh()[patchId_];
-    cellOwners_ = patch.faceCells();
+    patchInjectionBase::updateMesh(this->owner().mesh());
 }
 
 
@@ -155,6 +117,36 @@ Foam::scalar Foam::PatchFlowRateInjection<CloudType>::timeEnd() const
 }
 
 
+template<class CloudType>
+Foam::scalar Foam::PatchFlowRateInjection<CloudType>::flowRate() const
+{
+   const polyMesh& mesh = this->owner().mesh();
+
+    const surfaceScalarField& phi =
+        mesh.lookupObject<surfaceScalarField>(phiName_);
+
+    const scalarField& phip = phi.boundaryField()[patchId_];
+
+    scalar flowRateIn = 0.0;
+    if (phi.dimensions() == dimVelocity*dimArea)
+    {
+        flowRateIn = max(0.0, -sum(phip));
+    }
+    else
+    {
+        const volScalarField& rho =
+            mesh.lookupObject<volScalarField>(rhoName_);
+        const scalarField& rhop = rho.boundaryField()[patchId_];
+
+        flowRateIn = max(0.0, -sum(phip/rhop));
+    }
+
+    reduce(flowRateIn, sumOp<scalar>());
+
+    return flowRateIn;
+}
+
+
 template<class CloudType>
 Foam::label Foam::PatchFlowRateInjection<CloudType>::parcelsToInject
 (
@@ -164,10 +156,11 @@ Foam::label Foam::PatchFlowRateInjection<CloudType>::parcelsToInject
 {
     if ((time0 >= 0.0) && (time0 < duration_))
     {
-        scalar nParcels = fraction_*(time1 - time0)*parcelsPerSecond_;
+        scalar dt = time1 - time0;
 
-        cachedRandom& rnd = this->owner().rndGen();
+        scalar c = concentration_.value(0.5*(time0 + time1));
 
+        scalar nParcels = fraction_*parcelConcentration_*c*flowRate()*dt;
         label nParcelsToInject = floor(nParcels);
 
         // Inject an additional parcel with a probability based on the
@@ -177,7 +170,7 @@ Foam::label Foam::PatchFlowRateInjection<CloudType>::parcelsToInject
             nParcelsToInject > 0
          && (
                nParcels - scalar(nParcelsToInject)
-             > rnd.position(scalar(0), scalar(1))
+             > this->owner().rndGen().position(scalar(0), scalar(1))
             )
         )
         {
@@ -204,38 +197,11 @@ Foam::scalar Foam::PatchFlowRateInjection<CloudType>::volumeToInject
 
     if ((time0 >= 0.0) && (time0 < duration_))
     {
-        const polyMesh& mesh = this->owner().mesh();
+        scalar c = concentration_.value(0.5*(time0 + time1));
 
-        const surfaceScalarField& phi =
-            mesh.lookupObject<surfaceScalarField>(phiName_);
-
-        const scalarField& phip = phi.boundaryField()[patchId_];
-
-        scalar carrierVolume = 0.0;
-        if (phi.dimensions() == dimVelocity*dimArea)
-        {
-            const scalar flowRateIn = max(0.0, -sum(phip));
-            U0_ = -patchNormal_*flowRateIn/patchArea_;
-            carrierVolume = (time1 - time0)*flowRateIn;
-        }
-        else
-        {
-            const volScalarField& rho =
-                mesh.lookupObject<volScalarField>(rhoName_);
-            const scalarField& rhop = rho.boundaryField()[patchId_];
-
-            const scalar flowRateIn = max(0.0, -sum(phip/rhop));
-            U0_ = -patchNormal_*flowRateIn/patchArea_;
-            carrierVolume = (time1 - time0)*flowRateIn;
-        }
-
-        const scalar newParticles = concentration_*carrierVolume;
-
-        volume = pMeanVolume_*newParticles;
+        volume = c*(time1 - time0)*flowRate();
     }
 
-    reduce(volume, sumOp<scalar>());
-
     this->volumeTotal_ = volume;
     this->massTotal_ = volume*this->owner().constProps().rho0();
 
@@ -255,39 +221,15 @@ void Foam::PatchFlowRateInjection<CloudType>::setPositionAndCell
     label& tetPtI
 )
 {
-    if (cellOwners_.size() > 0)
-    {
-        cachedRandom& rnd = this->owner().rndGen();
-
-        label cellI = rnd.position<label>(0, cellOwners_.size() - 1);
-
-        cellOwner = cellOwners_[cellI];
-
-        // The position is between the face and cell centre, which could be
-        // in any tet of the decomposed cell, so arbitrarily choose the first
-        // face of the cell as the tetFace and the first point after the base
-        // point on the face as the tetPt.  The tracking will pick the cell
-        // consistent with the motion in the firsttracking step.
-        tetFaceI = this->owner().mesh().cells()[cellOwner][0];
-        tetPtI = 1;
-
-        // position perturbed between cell and patch face centres
-        const vector& pc = this->owner().mesh().C()[cellOwner];
-        const vector& pf =
-            this->owner().mesh().Cf().boundaryField()[patchId_][cellI];
-
-        const scalar a = rnd.sample01<scalar>();
-        const vector d = pf - pc;
-        position = pc + 0.5*a*d;
-    }
-    else
-    {
-        cellOwner = -1;
-        tetFaceI = -1;
-        tetPtI = -1;
-        // dummy position
-        position = pTraits<vector>::max;
-    }
+    patchInjectionBase::setPositionAndCell
+    (
+        this->owner().mesh(),
+        this->owner().rndGen(),
+        position,
+        cellOwner,
+        tetFaceI,
+        tetPtI
+    );
 }
 
 
@@ -300,8 +242,8 @@ void Foam::PatchFlowRateInjection<CloudType>::setProperties
     typename CloudType::parcelType& parcel
 )
 {
-    // set particle velocity
-    parcel.U() = U0_;
+    // set particle velocity to carrier velocity
+    parcel.U() = this->owner().U()[parcel.cell()];
 
     // set particle diameter
     parcel.d() = sizeDistribution_->sample();
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchFlowRateInjection/PatchFlowRateInjection.H b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchFlowRateInjection/PatchFlowRateInjection.H
index e6822a76d7..c666a92648 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchFlowRateInjection/PatchFlowRateInjection.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchFlowRateInjection/PatchFlowRateInjection.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -32,10 +32,10 @@ Description
       - Total mass to inject
       - Name of patch
       - Injection duration
-      - Initial parcel velocity
       - Injection target concentration/carrier volume flow rate
+    - Initial parcel velocity given by local flow velocity
     - Parcel diameters obtained by distribution model
-    - Parcels injected at cell centres adjacent to patch
+    - Parcels injected randomly across the patch
 
 SourceFiles
     PatchFlowRateInjection.C
@@ -46,15 +46,14 @@ SourceFiles
 #define PatchFlowRateInjection_H
 
 #include "InjectionModel.H"
+#include "patchInjectionBase.H"
+#include "TimeDataEntry.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-template<class Type>
-class TimeDataEntry;
-
 class distributionModel;
 
 /*---------------------------------------------------------------------------*\
@@ -64,22 +63,11 @@ class distributionModel;
 template<class CloudType>
 class PatchFlowRateInjection
 :
-    public InjectionModel<CloudType>
+    public InjectionModel<CloudType>,
+    public patchInjectionBase
 {
     // Private data
 
-        //- Name of patch
-        const word patchName_;
-
-        //- Id of patch
-        const label patchId_;
-
-        //- Patch area
-        scalar patchArea_;
-
-        //- Patch normal direction
-        vector patchNormal_;
-
         //- Name of carrier (mass or volume) flux field
         const word phiName_;
 
@@ -89,27 +77,15 @@ class PatchFlowRateInjection
         //- Injection duration [s]
         scalar duration_;
 
-        //- Concentration of particles to carrier [] (particles/m3)
-        const scalar concentration_;
+        //- Concentration profile of particle volume to carrier volume [-]
+        const TimeDataEntry<scalar> concentration_;
 
-        //- Number of parcels to introduce per second []
-        const label parcelsPerSecond_;
-
-        //- Initial parcel velocity [m/s]
-        vector U0_;
+        //- Parcels to introduce per unit volume flow rate m3 [n/m3]
+        const scalar parcelConcentration_;
 
         //- Parcel size distribution model
         const autoPtr<distributionModels::distributionModel> sizeDistribution_;
 
-        //- List of cell labels corresponding to injector positions
-        labelList cellOwners_;
-
-        //- Fraction of injection controlled by this processor
-        scalar fraction_;
-
-        //- Mean particle volume TODO: temporary measure - return from PDF
-        scalar pMeanVolume_;
-
 
 public:
 
@@ -152,6 +128,9 @@ public:
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Return the total volumetric flow rate across the patch [m3/s]
+        virtual scalar flowRate() const;
+
         //- Number of parcels to introduce relative to SOI
         virtual label parcelsToInject(const scalar time0, const scalar time1);