InflationInjection.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2011-2017 OpenFOAM Foundation
    Copyright (C) 2015-2022 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 "InflationInjection.H"
#include "mathematicalConstants.H"
#include "bitSet.H"
#include "cellSet.H"
#include "ListListOps.H"
 
using namespace Foam::constant::mathematical;
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
template<class CloudType>
Foam::InflationInjection<CloudType>::InflationInjection
(
    const dictionary& dict,
    CloudType& owner,
    const word& modelName
)
:
    InjectionModel<CloudType>(dict, owner, modelName, typeName),
    generationSetName_(this->coeffDict().lookup("generationCellSet")),
    inflationSetName_(this->coeffDict().lookup("inflationCellSet")),
    generationCells_(),
    inflationCells_(),
    duration_(this->coeffDict().getScalar("duration")),
    flowRateProfile_
    (
        Function1<scalar>::New
        (
            "flowRateProfile",
            this->coeffDict(),
            &owner.mesh()
        )
    ),
    growthRate_
    (
        Function1<scalar>::New
        (
            "growthRate",
            this->coeffDict(),
            &owner.mesh()
        )
    ),
    newParticles_(),
    volumeAccumulator_(0.0),
    fraction_(1.0),
    selfSeed_(this->coeffDict().getOrDefault("selfSeed", false)),
    dSeed_(SMALL),
    sizeDistribution_
    (
        distributionModel::New
        (
            this->coeffDict().subDict("sizeDistribution"),
            owner.rndGen()
        )
    )
{
    // Convert from user time to reduce the number of time conversion calls
    const Time& time = owner.db().time();
    duration_ = time.userTimeToTime(duration_);
    flowRateProfile_->userTimeToTime(time);
    growthRate_->userTimeToTime(time);
 
    if (selfSeed_)
    {
        this->coeffDict().readEntry("dSeed", dSeed_);
    }
 
    cellSet generationCells(this->owner().mesh(), generationSetName_);
 
    generationCells_ = generationCells.toc();
 
    cellSet inflationCells(this->owner().mesh(), inflationSetName_);
 
    // Union of cellSets
    inflationCells |= generationCells;
 
    inflationCells_ = inflationCells.toc();
 
    if (Pstream::parRun())
    {
        scalar generationVolume = 0.0;
 
        forAll(generationCells_, gCI)
        {
            label cI = generationCells_[gCI];
 
            generationVolume += this->owner().mesh().cellVolumes()[cI];
        }
 
        scalar totalGenerationVolume = generationVolume;
 
        reduce(totalGenerationVolume, sumOp<scalar>());
 
        fraction_ = generationVolume/totalGenerationVolume;
    }
 
    // Set total volume/mass to inject
    this->volumeTotal_ = fraction_*flowRateProfile_->integrate(0.0, duration_);
    this->massTotal_ *= fraction_;
}
 
 
template<class CloudType>
Foam::InflationInjection<CloudType>::InflationInjection
(
    const Foam::InflationInjection<CloudType>& im
)
:
    InjectionModel<CloudType>(im),
    generationSetName_(im.generationSetName_),
    inflationSetName_(im.inflationSetName_),
    generationCells_(im.generationCells_),
    inflationCells_(im.inflationCells_),
    duration_(im.duration_),
    flowRateProfile_(im.flowRateProfile_.clone()),
    growthRate_(im.growthRate_.clone()),
    newParticles_(im.newParticles_),
    volumeAccumulator_(im.volumeAccumulator_),
    fraction_(im.fraction_),
    selfSeed_(im.selfSeed_),
    dSeed_(im.dSeed_),
    sizeDistribution_(im.sizeDistribution_.clone())
{}
 
 
// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
template<class CloudType>
Foam::InflationInjection<CloudType>::~InflationInjection()
{}
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
template<class CloudType>
void Foam::InflationInjection<CloudType>::updateMesh()
{}
 
 
template<class CloudType>
Foam::scalar Foam::InflationInjection<CloudType>::timeEnd() const
{
    return this->SOI_ + duration_;
}
 
 
template<class CloudType>
Foam::label Foam::InflationInjection<CloudType>::parcelsToInject
(
    const scalar time0,
    const scalar time1
)
{
    const polyMesh& mesh = this->owner().mesh();
 
    List<DynamicList<typename CloudType::parcelType*>>& cellOccupancy =
        this->owner().cellOccupancy();
 
    scalar gR = growthRate_->value(time1);
 
    scalar dT = time1 - time0;
 
    // Inflate existing particles
 
    forAll(inflationCells_, iCI)
    {
        label cI = inflationCells_[iCI];
 
        typename CloudType::parcelType* pPtr = nullptr;
 
        forAll(cellOccupancy[cI], cPI)
        {
            pPtr = cellOccupancy[cI][cPI];
 
            scalar dTarget = pPtr->dTarget();
 
            pPtr->d() = min(dTarget, pPtr->d() + gR*dT);
        }
    }
 
    // Generate new particles
 
    newParticles_.clear();
 
    Random& rnd = this->owner().rndGen();
 
    // Diameter factor, when splitting particles into 4, this is the
    // factor that modifies the diameter.
    scalar dFact = sqrt(2.0)/(sqrt(3.0) + sqrt(2.0));
 
    if ((time0 >= 0.0) && (time0 < duration_))
    {
        volumeAccumulator_ +=
            fraction_*flowRateProfile_->integrate(time0, time1);
    }
 
    labelHashSet cellCentresUsed;
 
    // Loop escape counter
    label maxIterations = max
    (
        1,
        (10*volumeAccumulator_)
       /CloudType::parcelType::volume(sizeDistribution_().minValue())
    );
 
    label iterationNo = 0;
 
    // Info<< "Accumulated volume to inject: "
    //     << returnReduce(volumeAccumulator_, sumOp<scalar>()) << endl;
 
    while (!generationCells_.empty() && volumeAccumulator_ > 0)
    {
        if (iterationNo > maxIterations)
        {
            WarningInFunction
                << "Maximum particle split iterations ("
                << maxIterations << ") exceeded" << endl;
 
            break;
        }
 
        label cI = generationCells_
        [
            rnd.position
            (
                label(0),
                generationCells_.size() - 1
            )
        ];
 
        // Pick a particle at random from the cell - if there are
        // none, insert one at the cell centre.  Otherwise, split an
        // existing particle into four new ones.
 
        if (cellOccupancy[cI].empty())
        {
            if (selfSeed_ && !cellCentresUsed.found(cI))
            {
                scalar dNew = sizeDistribution_().sample();
 
                newParticles_.append
                (
                    vectorPairScalarPair
                    (
                        Pair<vector>(mesh.cellCentres()[cI], Zero),
                        Pair<scalar>(dSeed_, dNew)
                    )
                );
 
                volumeAccumulator_ -= CloudType::parcelType::volume(dNew);
 
                cellCentresUsed.insert(cI);
            }
        }
        else
        {
            label cPI = rnd.position(label(0), cellOccupancy[cI].size() - 1);
 
            // This has to be a reference to the pointer so that it
            // can be set to nullptr when the particle is deleted.
            typename CloudType::parcelType*& pPtr = cellOccupancy[cI][cPI];
 
            if (pPtr != nullptr)
            {
                scalar pD = pPtr->d();
 
                // Select bigger particles by preference
                if ((pD/pPtr->dTarget()) < rnd.sample01<scalar>())
                {
                    continue;
                }
 
                const point pP(pPtr->position());
                const vector& pU = pPtr->U();
 
                // Generate a tetrahedron of new positions with the
                // four new spheres fitting inside the old one, where
                // a is the diameter of the new spheres, and is
                // related to the diameter of the enclosing sphere, A,
                // by a = sqrt(2)*A/(sqrt(3) + sqrt(2));
 
                // Positions around the origin, which is the
                // tetrahedron centroid (centre of old sphere).
 
                // x = a/sqrt(3)
                // r = a/(2*sqrt(6))
                // R = sqrt(3)*a/(2*sqrt(2))
                // d = a/(2*sqrt(3))
 
                // p0(x, 0, -r)
                // p1(-d, a/2, -r)
                // p2(-d, -a/2, -r)
                // p3(0, 0, R)
 
                scalar a = pD*dFact;
 
                scalar x = a/sqrt(3.0);
                scalar r = a/(2.0*sqrt(6.0));
                scalar R = sqrt(3.0)*a/(2.0*sqrt(2.0));
                scalar d = a/(2.0*sqrt(3.0));
 
                scalar dNew = sizeDistribution_().sample();
                scalar volNew = CloudType::parcelType::volume(dNew);
 
                newParticles_.append
                (
                    vectorPairScalarPair
                    (
                        Pair<vector>(vector(x, 0, -r)  + pP, pU),
                        Pair<scalar>(a, dNew)
                    )
                );
                volumeAccumulator_ -= volNew;
 
                dNew = sizeDistribution_().sample();
                newParticles_.append
                (
                    vectorPairScalarPair
                    (
                        Pair<vector>(vector(-d, a/2, -r) + pP, pU),
                        Pair<scalar>(a, dNew)
                    )
                );
                volumeAccumulator_ -= volNew;
 
                dNew = sizeDistribution_().sample();
                newParticles_.append
                (
                    vectorPairScalarPair
                    (
                        Pair<vector>(vector(-d, -a/2, -r) + pP, pU),
                        Pair<scalar>(a, dNew)
                    )
                );
                volumeAccumulator_ -= volNew;
 
                dNew = sizeDistribution_().sample();
                newParticles_.append
                (
                    vectorPairScalarPair
                    (
                        Pair<vector>(vector(0, 0, R) + pP, pU),
                        Pair<scalar>(a, dNew)
                    )
                );
                volumeAccumulator_ -= volNew;
 
                // Account for the lost volume of the particle which
                // is to be deleted
                volumeAccumulator_ += CloudType::parcelType::volume
                (
                    pPtr->dTarget()
                );
 
                this->owner().deleteParticle(*pPtr);
 
                pPtr = nullptr;
            }
        }
 
        iterationNo++;
    }
 
    if (Pstream::parRun())
    {
        List<List<vectorPairScalarPair>> gatheredNewParticles
        (
            Pstream::nProcs()
        );
 
        gatheredNewParticles[Pstream::myProcNo()] = newParticles_;
 
        // Gather data onto master
        Pstream::gatherList(gatheredNewParticles);
 
        // Combine
        List<vectorPairScalarPair> combinedNewParticles
        (
            ListListOps::combine<List<vectorPairScalarPair>>
            (
                gatheredNewParticles,
                accessOp<List<vectorPairScalarPair>>()
            )
        );
 
        if (Pstream::master())
        {
            newParticles_ = combinedNewParticles;
        }
 
        Pstream::broadcast(newParticles_);
    }
 
    return newParticles_.size();
}
 
 
template<class CloudType>
Foam::scalar Foam::InflationInjection<CloudType>::volumeToInject
(
    const scalar time0,
    const scalar time1
)
{
    if ((time0 >= 0.0) && (time0 < duration_))
    {
        return fraction_*flowRateProfile_->integrate(time0, time1);
    }
 
    return 0.0;
}
 
 
template<class CloudType>
void Foam::InflationInjection<CloudType>::setPositionAndCell
(
    const label parcelI,
    const label,
    const scalar,
    vector& position,
    label& cellOwner,
    label& tetFacei,
    label& tetPti
)
{
    position = newParticles_[parcelI].first().first();
 
    this->findCellAtPosition
    (
        cellOwner,
        tetFacei,
        tetPti,
        position,
        false
    );
}
 
 
template<class CloudType>
void Foam::InflationInjection<CloudType>::setProperties
(
    const label parcelI,
    const label,
    const scalar,
    typename CloudType::parcelType& parcel
)
{
    parcel.U() = newParticles_[parcelI].first().second();
 
    parcel.d() = newParticles_[parcelI].second().first();
 
    parcel.dTarget() = newParticles_[parcelI].second().second();
}
 
 
template<class CloudType>
bool Foam::InflationInjection<CloudType>::fullyDescribed() const
{
    return false;
}
 
 
template<class CloudType>
bool Foam::InflationInjection<CloudType>::validInjection(const label)
{
    return true;
}
 
 
// ************************************************************************* //