diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectedParticleDistributionInjection/InjectedParticleDistributionInjection.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectedParticleDistributionInjection/InjectedParticleDistributionInjection.C new file mode 100644 index 0000000000000000000000000000000000000000..8254aec7ae18dba59683b0ca9d65147f3d31a2e6 --- /dev/null +++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectedParticleDistributionInjection/InjectedParticleDistributionInjection.C @@ -0,0 +1,525 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd. + \\/ 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 "InjectedParticleDistributionInjection.H" +#include "mathematicalConstants.H" +#include "PackedBoolList.H" +#include "injectedParticleCloud.H" + +using namespace Foam::constant; + +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + +template<class CloudType> +void Foam::InjectedParticleDistributionInjection<CloudType>::initialise() +{ + injectedParticleCloud ipCloud(this->owner().mesh(), cloudName_); + + List<label> tag(ipCloud.size()); + List<point> position(ipCloud.size()); + List<vector> U(ipCloud.size()); + List<scalar> soi(ipCloud.size()); + List<scalar> d(ipCloud.size()); + + // Flatten all data + label particlei = 0; + forAllConstIter(injectedParticleCloud, ipCloud, iter) + { + const injectedParticle& p = iter(); + tag[particlei] = p.tag(); + position[particlei] = p.position(); + U[particlei] = p.U(); + soi[particlei] = p.soi(); + d[particlei] = p.d(); + particlei++; + } + + // Combine all proc data + if (Pstream::parRun()) + { + List<List<label>> procTag(Pstream::nProcs()); + procTag[Pstream::myProcNo()].transfer(tag); + Pstream::gatherList(procTag); + Pstream::scatterList(procTag); + tag = + ListListOps::combine<List<label>> + ( + procTag, accessOp<List<label>>() + ); + + List<List<point>> procPosition(Pstream::nProcs()); + procPosition[Pstream::myProcNo()].transfer(position); + Pstream::gatherList(procPosition); + Pstream::scatterList(procPosition); + position = + ListListOps::combine<List<point>> + ( + procPosition, accessOp<List<point>>() + ); + + List<List<vector>> procU(Pstream::nProcs()); + procU[Pstream::myProcNo()].transfer(U); + Pstream::gatherList(procU); + Pstream::scatterList(procU); + U = + ListListOps::combine<List<vector>> + ( + procU, accessOp<List<vector>>() + ); + + List<List<scalar>> procSOI(Pstream::nProcs()); + procSOI[Pstream::myProcNo()].transfer(soi); + Pstream::gatherList(procSOI); + Pstream::scatterList(procSOI); + soi = + ListListOps::combine<List<scalar>> + ( + procSOI, accessOp<List<scalar>>() + ); + + List<List<scalar>> procD(Pstream::nProcs()); + procD[Pstream::myProcNo()].transfer(d); + Pstream::gatherList(procD); + Pstream::scatterList(procD); + d = + ListListOps::combine<List<scalar>> + ( + procD, accessOp<List<scalar>>() + ); + } + + label maxTag = -1; + forAll(tag, particlei) + { + maxTag = max(maxTag, tag[particlei]); + } + + label nInjectors = maxTag + 1; + List<scalar> injStartTime(nInjectors, GREAT); + List<scalar> injEndTime(nInjectors, -GREAT); + List<DynamicList<point>> injPosition(nInjectors); + List<DynamicList<vector>> injU(nInjectors); + List<DynamicList<scalar>> injDiameter(nInjectors); + + // Cache the particle information per tag + forAll(tag, i) + { + const label tagi = tag[i]; + const scalar t = soi[i]; + injStartTime[tagi] = min(t, injStartTime[tagi]); + injEndTime[tagi] = max(t, injEndTime[tagi]); + injPosition[tagi].append(position[i]); + injU[tagi].append(U[i]); + injDiameter[tagi].append(d[i]); + } + + // Remove single particles and injectors where injection interval is 0 + // - cannot generate a volume flow rate + scalar sumVolume = 0; + startTime_.setSize(nInjectors, 0); + endTime_.setSize(nInjectors, 0); + sizeDistribution_.setSize(nInjectors); + position_.setSize(nInjectors); + U_.setSize(nInjectors); + volumeFlowRate_.setSize(nInjectors, 0); + + scalar minTime = GREAT; + + // Populate injector properties, filtering out invalid entries + cachedRandom& rnd = this->owner().rndGen(); + label injectori = 0; + forAll(injDiameter, i) + { + const DynamicList<scalar>& diameters = injDiameter[i]; + const label nParticle = diameters.size(); + const scalar dTime = injEndTime[i] - injStartTime[i]; + + if ((nParticle > 1) && (dTime > ROOTVSMALL)) + { + minTime = min(minTime, injStartTime[i]); + + startTime_[injectori] = injStartTime[i]; + endTime_[injectori] = injEndTime[i]; + + // Re-sample the cloud data + position_[injectori].setSize(resampleSize_); + U_[injectori].setSize(resampleSize_); + List<point>& positioni = position_[injectori]; + List<vector>& Ui = U_[injectori]; + + for (label samplei = 0; samplei < resampleSize_; ++samplei) + { + label posi = rnd.globalPosition<label>(0, nParticle - 1); + positioni[samplei] = injPosition[i][posi] + positionOffset_; + Ui[samplei] = injU[i][posi]; + } + + // Calculate the volume flow rate + scalar sumPow3 = 0; + forAll(diameters, particlei) + { + sumPow3 += pow3(diameters[particlei]); + } + + const scalar volume = sumPow3*mathematical::pi/16.0; + sumVolume += volume; + volumeFlowRate_[injectori] = volume/dTime; + + // Create the size distribution using the user-specified bin width + sizeDistribution_.set + ( + injectori, + new distributionModels::general + ( + diameters, + binWidth_, + this->owner().rndGen() + ) + ); + + injectori++; + } + } + + // Resize + startTime_.setSize(injectori); + endTime_.setSize(injectori); + position_.setSize(injectori); + U_.setSize(injectori); + volumeFlowRate_.setSize(injectori); + sizeDistribution_.setSize(injectori); + + // Reset start time to zero + forAll(startTime_, injectori) + { + startTime_[injectori] -= minTime; + endTime_[injectori] -= minTime; + } + + + // Set the volume of parcels to inject + this->volumeTotal_ = sumVolume; + + // Provide some feedback + Info<< " Read " << position_.size() << " injectors with " + << tag.size() << " total particles" << endl; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template<class CloudType> +Foam::InjectedParticleDistributionInjection<CloudType>:: +InjectedParticleDistributionInjection +( + const dictionary& dict, + CloudType& owner, + const word& modelName +) +: + InjectionModel<CloudType>(dict, owner, modelName, typeName), + cloudName_(this->coeffDict().lookup("cloud")), + startTime_(this->template getModelProperty<scalarList>("startTime")), + endTime_(this->template getModelProperty<scalarList>("endTime")), + position_(this->template getModelProperty<List<vectorList>>("position")), + positionOffset_(this->coeffDict().lookup("positionOffset")), + volumeFlowRate_ + ( + this->template getModelProperty<scalarList>("volumeFlowRate") + ), + U_(this->template getModelProperty<List<vectorList>>("U")), + binWidth_(readScalar(this->coeffDict().lookup("binWidth"))), + sizeDistribution_(), + parcelsPerInjector_ + ( + ceil(readScalar(this->coeffDict().lookup("parcelsPerInjector"))) + ), + resampleSize_ + ( + this->coeffDict().template lookupOrDefault<label>("resampleSize", 100) + ), + applyDistributionMassTotal_ + ( + readBool(dict.lookup("applyDistributionMassTotal")) + ), + ignoreOutOfBounds_ + ( + this->coeffDict().lookupOrDefault("ignoreOutOfBounds", false) + ), + nParcelsInjected_(this->parcelsAddedTotal()), + nParcelsInjected0_(0), + currentInjectori_(0), + currentSamplei_(0) +{ + if (startTime_.size()) + { + // Restart + sizeDistribution_.setSize(startTime_.size()); + forAll(sizeDistribution_, i) + { + const word dictName("distribution" + Foam::name(i)); + dictionary dict; + this->getModelDict(dictName, dict); + + sizeDistribution_.set + ( + i, + new distributionModels::general(dict, this->owner().rndGen()) + ); + } + } + else + { + // Clean start + initialise(); + } + + // Set the mass of parcels to inject from distribution if requested + if (applyDistributionMassTotal_) + { + this->massTotal_ = this->volumeTotal_*this->owner().constProps().rho0(); + Info<< " Set mass to inject from distribution: " + << this->massTotal_ << endl; + } +} + + +template<class CloudType> +Foam::InjectedParticleDistributionInjection<CloudType>:: +InjectedParticleDistributionInjection +( + const InjectedParticleDistributionInjection<CloudType>& im +) +: + InjectionModel<CloudType>(im), + cloudName_(im.cloudName_), + startTime_(im.startTime_), + endTime_(im.endTime_), + position_(im.position_), + positionOffset_(im.positionOffset_), + volumeFlowRate_(im.volumeFlowRate_), + U_(im.U_), + binWidth_(im.binWidth_), + sizeDistribution_(im.sizeDistribution_.size()), + parcelsPerInjector_(im.parcelsPerInjector_), + resampleSize_(im.resampleSize_), + applyDistributionMassTotal_(im.applyDistributionMassTotal_), + ignoreOutOfBounds_(im.ignoreOutOfBounds_), + nParcelsInjected_(im.nParcelsInjected_), + nParcelsInjected0_(im.nParcelsInjected0_), + currentInjectori_(0), + currentSamplei_(0) +{ + forAll(sizeDistribution_, injectori) + { + if (sizeDistribution_.set(injectori)) + { + sizeDistribution_.set + ( + injectori, + new distributionModels::general + ( + im.sizeDistribution_[injectori] + ) + ); + } + } +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +template<class CloudType> +Foam::InjectedParticleDistributionInjection<CloudType>:: +~InjectedParticleDistributionInjection() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template<class CloudType> +void Foam::InjectedParticleDistributionInjection<CloudType>::updateMesh() +{} + + +template<class CloudType> +Foam::scalar +Foam::InjectedParticleDistributionInjection<CloudType>::timeEnd() const +{ + return max(endTime_); +} + + +template<class CloudType> +Foam::label +Foam::InjectedParticleDistributionInjection<CloudType>::parcelsToInject +( + const scalar time0, + const scalar time1 +) +{ + // Ensure all procs know the latest parcel count + nParcelsInjected_ += returnReduce(nParcelsInjected0_, sumOp<label>()); + nParcelsInjected0_ = 0; + + scalar targetVolume = 0; + forAll(startTime_, injectori) + { + if (time1 > startTime_[injectori]) + { + scalar totalDuration = + min(time1, endTime_[injectori]) - startTime_[injectori]; + + targetVolume += volumeFlowRate_[injectori]*totalDuration; + } + } + + const label targetParcels = + round + ( + scalar(startTime_.size()*parcelsPerInjector_) + *targetVolume/this->volumeTotal_ + ); + + const label nParcels = targetParcels - nParcelsInjected_; + + return nParcels; +} + + +template<class CloudType> +Foam::scalar +Foam::InjectedParticleDistributionInjection<CloudType>::volumeToInject +( + const scalar time0, + const scalar time1 +) +{ + scalar volume = 0; + forAll(startTime_, injectori) + { + if ((time1 > startTime_[injectori]) && (time1 <= endTime_[injectori])) + { + scalar duration = min(time1, endTime_[injectori]) - time0; + volume += volumeFlowRate_[injectori]*duration; + } + } + + return volume; +} + + +template<class CloudType> +void Foam::InjectedParticleDistributionInjection<CloudType>::setPositionAndCell +( + const label parcelI, + const label nParcels, + const scalar time, + vector& position, + label& cellOwner, + label& tetFaceI, + label& tetPtI +) +{ + cachedRandom& rnd = this->owner().rndGen(); + currentInjectori_ = rnd.globalPosition<label>(0, position_.size() - 1); + currentSamplei_ = rnd.globalPosition<label>(0, resampleSize_ - 1); + + position = position_[currentInjectori_][currentSamplei_]; + + // Cache all mesh props for each position? + this->findCellAtPosition + ( + cellOwner, + tetFaceI, + tetPtI, + position + ); +} + + +template<class CloudType> +void Foam::InjectedParticleDistributionInjection<CloudType>::setProperties +( + const label parcelI, + const label, + const scalar, + typename CloudType::parcelType& parcel +) +{ + // Set particle velocity + parcel.U() = U_[currentInjectori_][currentSamplei_]; + + // Set particle diameter + parcel.d() = sizeDistribution_[currentInjectori_].sample(); + + // Increment number of particles injected + // Note: local processor only! + nParcelsInjected0_++; +} + + +template<class CloudType> +bool +Foam::InjectedParticleDistributionInjection<CloudType>::fullyDescribed() const +{ + return false; +} + + +template<class CloudType> +bool Foam::InjectedParticleDistributionInjection<CloudType>::validInjection +( + const label +) +{ + return true; +} + + +template<class CloudType> +void Foam::InjectedParticleDistributionInjection<CloudType>::info(Ostream& os) +{ + InjectionModel<CloudType>::info(os); + + if (this->writeTime()) + { + this->setModelProperty("startTime", startTime_); + this->setModelProperty("endTime", endTime_); + this->setModelProperty("position", position_); + this->setModelProperty("volumeFlowRate", volumeFlowRate_); + this->setModelProperty("U", U_); + forAll(sizeDistribution_, i) + { + const distributionModels::general& dist = sizeDistribution_[i]; + const word dictName("distribution" + Foam::name(i)); + dictionary dict(dist.writeDict(dictName)); + this->setModelProperty(dictName, dict); + } + } +} + + +// ************************************************************************* // diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectedParticleDistributionInjection/InjectedParticleDistributionInjection.H b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectedParticleDistributionInjection/InjectedParticleDistributionInjection.H new file mode 100644 index 0000000000000000000000000000000000000000..3ccfa5db278227114c5ed3381a7ab980709eeb07 --- /dev/null +++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectedParticleDistributionInjection/InjectedParticleDistributionInjection.H @@ -0,0 +1,263 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2015-2016 OpenCFD Ltd. + \\/ 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::InjectedParticleDistributionInjection + +Description + Interrogates an injectedParticleCloud to convert the raw particle data + into a set of 'binned' injectors. + + The bins are set according to the particle \c tag property, from which: + - diameters are converted into \c general distributions with a + user-specified bin width + - raw velocity and diameter data are resampled and stored to provide + variations per injector + + The mass to inject can be set according to the raw input data mass total + by using the \c applyDistributionMassTotal switch + + +Usage + \verbatim + model1 + { + type injectedParticleDistributionInjection; + SOI 0; + parcelBasisType mass; + cloud eulerianParticleCloud; + positionOffset (-0.025 2 -0.025); + binWidth 0.1e-3; + parcelsPerInjector 500; + resampleSize 100; // optional + applyDistributionMassTotal yes; + + // Placeholder only when using applyDistributionMassTotal + massTotal 0; + } + \endverbatim + +Note + The each injector location is assumed to be operating under steady + conditions, i.e. using a constant flow rate profile + +SourceFiles + InjectedParticleDistributionInjection.C + +See also + Foam::injectedParticle + Foam::injectedParticleCloud + Foam::functionObjects::extractEulerianParticles + Foam::distributionModels::general + +\*---------------------------------------------------------------------------*/ + +#ifndef InjectedParticleDistributionInjection_H +#define InjectedParticleDistributionInjection_H + +#include "InjectionModel.H" +#include "general.H" +#include "Switch.H" +#include "vectorList.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class InjectedParticleDistributionInjection Declaration +\*---------------------------------------------------------------------------*/ + +template<class CloudType> +class InjectedParticleDistributionInjection +: + public InjectionModel<CloudType> +{ +protected: + + // Protected data + + //- Name of cloud used to seed the new particles + const word cloudName_; + + //- List of start time per injector + scalarList startTime_; + + //- List of end time per injector + scalarList endTime_; + + //- List of position per injector + List<vectorList> position_; + + //- Position offset to apply to input positions + vector positionOffset_; + + //- List of volume flow rate per injector [m3/s] + scalarList volumeFlowRate_; + + //- List of parcel velocity per injector + List<vectorList> U_; + + //- Bin width when generating particle distributions + scalar binWidth_; + + //- List of size distribution model per injector + PtrList<distributionModels::general> sizeDistribution_; + + //- Target number of parcels to inject per injector + scalar parcelsPerInjector_; + + //- Resample size + label resampleSize_; + + //- Flag to apply mass calculated from distribution instead of + // InjectionModel massTotal_ + bool applyDistributionMassTotal_; + + //- Flag to suppress errors if particle injection site is out-of-bounds + Switch ignoreOutOfBounds_; + + //- Running total of number of parcels injected + label nParcelsInjected_; + + //- Number of parcels injected in last step (local proc only) + label nParcelsInjected0_; + + //- Current injector + label currentInjectori_; + + //- Current sample + label currentSamplei_; + + + // Protected Member Functions + + //- Initialise injectors + void initialise(); + + +public: + + //- Runtime type information + TypeName("injectedParticleDistributionInjection"); + + + // Constructors + + //- Construct from dictionary + InjectedParticleDistributionInjection + ( + const dictionary& dict, + CloudType& owner, + const word& modelName + ); + + //- Construct copy + InjectedParticleDistributionInjection + ( + const InjectedParticleDistributionInjection<CloudType>& im + ); + + //- Construct and return a clone + virtual autoPtr<InjectionModel<CloudType> > clone() const + { + return autoPtr<InjectionModel<CloudType> > + ( + new InjectedParticleDistributionInjection<CloudType>(*this) + ); + } + + + //- Destructor + virtual ~InjectedParticleDistributionInjection(); + + + // Member Functions + + //- Set injector locations when mesh is updated + virtual void updateMesh(); + + //- Return the end-of-injection time + scalar timeEnd() const; + + //- Number of parcels to introduce relative to SOI + virtual label parcelsToInject(const scalar time0, const scalar time1); + + //- Volume of parcels to introduce relative to SOI + virtual scalar volumeToInject(const scalar time0, const scalar time1); + + + // Injection geometry + + //- Set the injection position and owner cell, tetFace and tetPt + virtual void setPositionAndCell + ( + const label parcelI, + const label nParcels, + const scalar time, + vector& position, + label& cellOwner, + label& tetFaceI, + label& tetPtI + ); + + //- Set the parcel properties + virtual void setProperties + ( + const label parcelI, + const label nParcels, + const scalar time, + typename CloudType::parcelType& parcel + ); + + //- Flag to identify whether model fully describes the parcel + virtual bool fullyDescribed() const; + + //- Return flag to identify whether or not injection of parcelI is + // permitted + virtual bool validInjection(const label parcelI); + + + // I-O + + //- Write injection info to stream + void info(Ostream& os); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "InjectedParticleDistributionInjection.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* //